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ABSTRACT 



A modified binary search routine was developed by 
Tysver to provide inputs (estimates of parameters of thres- 
hold stimulus distributions) to the Probit and Staircase 
techniques for sensitivity testing. A high speed digital 
computer was used to provide simulated sensitivity data 
in order to test the modified search (more specifically, 
study the parameter estimates generated by the search). 

The results show that the modified search does not 
give practical input information to the techniques mentioned 
above. However, further refinements of extensions to the 
basic search patterns should yield better estimates of the 
distribution parameters. 
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I . INTRODUCTION 



The estimations of parameters describing distributions 
of outcomes in various physical situations are usually 
easy to compute and have known accuracy. The best estimate 
of the mean time to failure of a light bulb, for example, 
is simply the average failure time of a sample of light 
bulbs. The best estimates of the mean and variance of a 
normally distributed outcome are the sample average and 
sample variance using well known formulas . The samples 
in these cases are made up of a collection of point esti- 
mates of the population mean. 

However, there exists a class of physical situations 
in nature for which the parameters of the distributions 
cannot be estimated by the usual point estimate techniques. 
Examples are the threshold stimulus to detonate an explosive 
charge, the amount of insecticide necessary to kill a pesky 
mosquito, and the intensity of a light source necessary 
for visual perception at a given range. These examples 
share a common characteristic: they all deal with an appli- 

cation of some form of impulse stimulus. The study of this 
class of physical situations has been labeled "sensitivity 
analysis . " 

When applying a stimulus to a certain sense, one of two 
responses can occur. The subject has a positive response 
(e.g., detonates, dies, sees the light source) or a 
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negative response (e.g., fails to detonate, lives, fails to 
see the light source). The response is simply a Bernoulli 
random variable and the probability of a response varies 
with the stimulus level applied. 

Major efforts in the past have been devoted to statis- 
tical techniques centered on the density function and point 
estimates. However, estimates of parameters of sensitivity 
distributions cannot be established by usual point estima- 
tion methods since data realized in sensitivity testing 
is obtained from distribution functions, not density 
functions . 

A few techniques have been developed over the years 
to accomplish the task of providing parameter estimates 
from sensitivity data. These include the Probit technique, 
as described by Finney £ 2 ] , the Staircase or Bruceton method 
reviewed by Dixon and Mood [lj , the Countback method of 
Lewis Ql] , and others. Many methods are devoted to finding 
the 50# level (the level at which one would expect 50 % 
positive responses to the stimulus [V) ) . Others claim to 
estimate parameters of sensitivity distributions accurately 
but require knowledge or assumptions about the parameters 
prior to testing. 

The Modified Binary Search (MBS) is a method designed 
by Tysver |jf] to examine physical situations for which sensi- 
tivity testing is applicable. Various levels of stimuli 
are used as a basis for computing estimates of parameters 
of the appropriate threshold stimulus distribution. 



6 



Tysver proposed that the estimates from MBS be used as 
inputs to previously mentioned tests (e.g. Probit, Stair- 
case) and also in predictions for safety and reliability. 

The results of were examined and tested through 

the use of simulation conducted on a high-speed electronic 
computer. Standard Monte Carlo techniques were used to 
generate the required sensitivity data. 
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II. THE MODIFIED BINARY SEARCH 



A. UNDERLYING SENSITIVITY MODEL 

A brief description of the model used in the Modified 
Binary Search (MBS) is as follows. Let x be an applied 
stimulus level (xe(0,°°)). Define Y = f(x) to be a Bernoulli 
random variable with realizations y = 1 for a positive 
response and y = 0 for no response. Next, define a response 
function p(x) where 

p(x) = Prob(Y = l|x), 

and assume that p(0) = 0 and p(°°) =1. An investigator, 
however, should be able to determine a shorter interval 
(a,b) with 0<a<x<b<°° where p(a) ~ 0 and p(b) = 1 
with a high degree of confidence. The response function 
is graphed in Figure 1. 



p(x) 




Figure 1. The Response Function 
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In other words, p(x) is the cumulative distribution func- 
tion for a random variable X (the threshold stimulus), 
where 

p(x) = Prob(X < x) . 

If the applied stimulus level is greater than the thres- 
hold stimulus, X, then the response will be y = 1. 

If x < X, then the response will be y = 0. Note that 

Prob(Y = l|x) = Prob(X < x) = p(x) 

and 

Prob(Y = 0|x) = Prob(X > x) = 1 - p(x). 

The threshold stimulus, X, or an appropriate transforma- 
tion of the same, is assumed to be normally distributed. 

The validity of this assumption does not appear to be 
critical in the vicinity of the 50% response level, especial- 
ly when the investigation is limited to small samples [6] . 

With the sensitivity model defined, the investigator 
is now faced with the task of preselecting the applied 
stimulus levels (xi, x 2 ,..., x n ) or determining a precedure 
to choose each stimulus level x^ +1 based upon previously 
tested stimulus levels (xi, x 2 ,..., x^) and the observed 
results (yi, y 2 ,..., y.^). Once a sampling procedure 
has been selected, he conducts the investigation. The 
sample data ( y i , y 2 , . . . , y ) he obtains contain information 
on the distribution function of the threshold stimulus. 

Each response (y^) tells the investigator whether or not 
the random threshold stimulus (X^) was greater than the 
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applied stimulus level (x 1 ). The investigator is thus 
faced with estimating parameters of the threshold stimulus' 
density function from data describing the distribution 
function . 

B. DESCRIPTION OF MBS 

The investigator's first task is to choose reasonable 

endpoints a, and b over the range of stimulus levels where 

p(a) 0 and p(b) ss 1. Sampling outside the interval 

(a, b) will be unnecessary if the endpoints are chosen 

properly. For his first trial he selects Xj = (a + b)/2. 

If y = 1, then he selects x 2 = (a + b)/4 for his next 

trial. If y = 0, however, he selects x 2 = 3(a + b)/4. 

The complete search procedure is diagramed in Figures 2 

through 5. If it appears that the search is converging 

to one of the endpoints, then the investigator should 

question his initial choice of that endpoint. 

As in other sensitivity analysis techniques, MBS does 

not terminate until an "inversion" occurs (y^ = 1 and 

y. = 0 for some i, j where x. < x.). Six equidistant 
J J- J 

levels are tested in the vicinity of the inversion. 

Let the lowest of the six levels be denoted by x^ and the 
highest by x . Define the interval (x T , x ) as the "pri- 
mary region" and the distance separating any two adjacent 
levels in this primary region as "DELX." Any levels tested 
outside of the primary region are separated by at least 
DELX. The preliminary search ends with one of three pos- 
sible sequential outcomes, Si, Sz, or S3. The number of 
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trials necessary to produce a sequence may vary with each 
experiment, but the primary region will consist of only 
three possible outcomes as depicted in Figure 6. All 
levels tested below the primary region yielded "0" respon- 
ses while those tested above the region yielded positive 
responses . 

y = 0 . . . .001011. . . .1 

51 • 

x = a. . . .xxxxxx. . . . b 

y = 0. . . .101011. . . .1 

52 t — ' ■ 

x = a. . . .xxxxxx. . . .b 

y = 0. . . .001010. . . .1 

53 : - — 

x = a. . . .xxxxxx. . . .b 

Figure 6. Primary Region 

Once a terminal sequence (Si , S 2 , or S 3 ) has been 
obtained, maximum likelihood estimates of the parameters 
of the stimulus distribution are obtained. 

Let n be the number of trials in a terminal sequence. 
The estimates, y and a, are found by maximizing: 

n 

n Prob^ = yjx^ b,o) Equation (1) 

i=l 

where 

f N(x 1 ; y , a ) if y ± = 1 

Prob(Y 1 = y i |x i ; y,a) = | 

I 1 - N^j y , a ) if y 1 = 0 
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and 



N(x; y,a) = Prob(X i x; y,c) = 



(x-y ) : 
2a 2 



/2?a 



dx . 



In other words, given a sequence of responses (yi , 
y 2 , . . . , y ) for stimulus levels (x x , x 2 , . . . , x n )> find 
the maximum likelihood estimates (mle's) of the parameters 
that describe the normal distribution from which the data 
was realized. These estimates, y and a, can then be used 
as inputs to the Staircase or Probit techniques. If more 
trials are desired in order to yield more accurate esti- 
mates, then the investigator should follow the extension 
procedures as described below. 



C. EXTENSION OP SEARCH 
1 . Extension of S’i 

For the terminal sequence Si the responses are 
fairly symmetric about y, so that additional trials should 
be carried out in the primary region. Figure 7 depicts 
the stimulus levels to test in the event three or five 
additional trials are desired. 

0 0 10 11 

xxx 

Three Trials 

0 0 10 11 

X X X X X 

Five Trials 

Figure 7. Extension of Si 
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2 . Extension of S2 or S3 



It appears that one of the ends of the primary 
region is fairly well pinned down while the other is not. 
Therefore additional trials are recommended as shown in 
Figure 8. 

10 10 11 

xxx 

5 2 Extension 

0 0 10 10 

xxx 

5 3 Extension 

Figure 8. Extensions of S 2 and S 3 

D. EXAMPLE OF A SEARCH USING MBS 

A sample experiment may go as follows. Suppose an 
investigator is testing for the flashpoint of an explosive. 
This type of an event certainly qualifies as belonging 
to the sensitivity class. From past experience he roughly 

o o 

chooses a = 0 C and b = 1000 C. His first trial is con- 

O 

ducted at 500 C and he observes that the explosive ignites. 
The next level (temperature) at which to test the explosive 

O 

is thus 250 C. Continuing his experiment he arrives at 
terminal sequence Si after eight trials. 



Trial 


Level ( °C ) 


Outcome 


1 


500 


1 (ignition) 


2 


250 


0 (no ignition) 


3 


750 


1 (s*) 


4 


375 


0 
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Trial 


Level ( °C) 


Outcome 


5 


625 


1 (s*) 


6 


437.5 


0 


7 


562.5 


0 (inversion) 


8 


687.5 


1 



An ordered listing of the trial levels and respective 
outcomes is as follows. The equations for the mle ’ s will 
appear later. 



y = 0 



x = 0 



250 375 



500 



625 



750 



437.5 562.5 687.5 



1 

1000 (°C) 



A o 

The maximum likelihood estimates are: u = 531.25 C; 
a = 81.25 °C. 
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III. NATURE OF PROBLEM 



In order to adequately test the MBS a digital computer 
was used for both simulation and data processing. It 
appeared that a standard analytical method (e.g. using 
derivatives) to locate the maximum of Equation (1) was 
not feasible. However, a favorable property of Equation (1) 
is that the function is unimodal over the entire range of 
y and a. Thus, search routines (e.g. Direct, Hookes-Jeeves ) 
can be used to find estimates of the maximum of Equation (1) 
and the associated y and a as well. 

All estimates of threshold parameters calculated in 
[5] were done by hand, or desk calculator, thus permitting 
the possibility of significant round off error. The only 
inputs (response data)- to Equation (1) were taken from the 
primary inversion region, while responses outside this 
region were ignored. This seemed to have the greatest 
effect of the sequences S 2 and S 3 . Inclusion of all res- 
ponses was considered reasonable and thus written into 
the computer program used in this study. 

In an earlier study on MBS Hicks [3] claimed that o 
underestimated a significantly with no predictable bias. 

The distribution of a was therefore investigated in this 
study for each of the terminal sequences. In addition, 
the distribution of y was investigated. 
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IV. SIMULATION 



A. DESCRIPTION 

The computer program was written in the Fortran IV 
programming language and executed on an IBM 360/67 computer. 
The program was written in two sections, I and II. 

Section I conducted a simulated sensitivity experiment 
leading to a terminal sequence (Sj , S 2 , or S 3 ) . A standard 
normal distribution was used as the source of the sensiti- 
vity data. The distribution endpoints (a, b) were varied 
in order to test MBS over a wide range of starting points. 

In most cases a~U(-9, -3) and b~U(3> 9). At each stimulus 
level, x^ , a random number was drawn from the simulated 
distribution and tested. If the random number was greater 
than x^ , then y^ = 0. The program followed the complete 
search procedure as illustrated in Figures 2 through 5* 

A simple flow chart is diagrammed in Figure 9* 

Section II took the vectors (x 1} x 2 ,..., x ) and 
(y t , y 2 , . . . , y ) from Section I and calculated estimates 
of y and a to within .03 DELX of the values y* and a* 
that maximize Equation (1). Initially a "direct" search 
was employed to locate the mle ' s for a given sequence. 
Different guesses of y and a were used as inputs to Equa- 
tion (1). This method proved to be very inefficient. 

A Hookes-Jeeves search routine for two variables was tried 
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with much better success (a copy of this routine can be 
found at the end of the thesis). 



B. RESULTS 

1 . Verification of Previous Estimates 

The majority of estimates in [ 5 ] were verified by 

the computer simulations. Table 1 lists the results for 

comparison. The values of y are listed in displacement 

from x„ where x is the smallest level at which a positive 
s s 

response occurred. This displacement is given in multiples 
of DELX (dx). The values of a are also given as multiples 
of DELX (dx). 

The estimates, y and a, for Si are not affected 
by the total number of levels tested. However, the esti- 
mates from an S 2 or S 3 do depend on the total number of 
trials in the sequence. The values for y and a for S 2 
and S3 listed in Table 1 are valid for sequences of eight 
or more trials. If an investigator reaches S 2 or S 3 with 
seven or fewer trials he should extend his search or repeat 
it altogether. 

The frequency of occurrences of Si, S 2 , and S 3 
from a sample of 2885 sequences was: 



51 

5 2 

5 3 



frequency 

2226 

330 

329 



percentage 

77.2 

11.4 

11.4 
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Revised Estimates 



Sequence Previous Estimates [ 5 ] 



A A A A 

y a y 0 



Si 


x + dx/2 
s 


1 . 3dx 


x + . 5dx 
s 


1 . 3dx 


S 2 


x - dx/4 
s 


6dx 


x +1.2dx 
s 


3dx 


S 3 


x +13dx/4 

5 


6dx 


x +1.8dx 
s 


3dx 






Si extension 


- 






(three 


additional trials) 






A 


A 


/V 


A 


outcome 


y 


a 


y 


0 


0 0 0 


x s +l . 5dx 


1.5dx 


x + . 9dx 
s 


1 . 5dx 


0 0 1 


x +.75dx 

0 


dx 


x g + . 2dx 


dx 


0 10 


V dx 


2dx 


x g + . 3dx 


1 . 6dx 


Oil 


x + . 25dx 
s 


dx 


x - . 2dx 
s 


dx 


10 0 


V dx 


2.75dx 


x s + . 3dx 


2 . 3dx 


10 1 


x s 


2dx 


x g - . 3dx 


1 . 6dx 


110 


x s 


2 . 75dx 


x s -.3dx 


2 . 3dx 


111 


x - dx 

s 


1 . 5dx 


x s -.9dx 


1 . 5dx 



dx E DELX 
TABLE 1 

SIMULATION RESULTS 
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2. Biased Estimates 



Levels tested in the vicinity of the mean of the 
stimulus distribution offer the most informative data, 
while those tested within a few standard deviations of 
the mean aid in pinning down the tails of the tlistribution . 
The Probit technique requires the investigation to com- 
mence fairly close to the mean. The Probit and Staircase 
methods use a constant step size between levels and thus 
perform best when the step size x c.' For instance, if the 
step size >>o, then an investigator might never test within 
three or four standard deviations of the mean. If the 
step size <<a, then one would need many trials to insure 
that an informative range of levels was tested in the 
vicinity of the mean. A useful estimate of a is one that 
falls within the range of a/2 to 2a. 

The values listed in Table 1 are functions of x 

s 

and DELX, or DELX alone. The level x (smallest level 
at which a positive response occurs) is random but tends 
toward the mean of the distribution. DELX is also random 
but depends upon the selection of the endpoints, a and b, 
the threshold variance, and the sequential outcomes 
(y i , y 2 , . . . , y ) . Thus DELX can take on values over a 
wide range implying y and $ can as well. 

On the basis of substantial sampling, estimates 
of y appear to be unbiased for Sj, but seem to have a 
small predictable bias for S 2 and S 3 . The average u from 
simulated N(0, 1) sensitivity data was about -0.1 for S 2 
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and +0.1 for S 3 . This would correspond to p - .la and y + .la 
for N (y,a) data. 

The estimate a is on the other hand not so well 
behaved. It appears to have unpredictable bias for all 
sequences (Hicks [3] ) . Less than half of the estimates 
of a generated from hundreds of simulated trials were within 
the desired range (a/2 to 2a). The majority were less 
than a/2. 

Of primary concern is why the apparent bias? Is 
Equation (1) in error? Equation (1) does not depend upon 
the number of trials necessary for an inversion, whereas 
DELX does. It appears that MBS generates too many trials 
(i.e. DELX is frequently small in comparison to a when the 
first inversion occurs). Hence, the average estimate of 
a is much smaller than the true standard deviation. 

Consider an estimate a generated from an Sj sequence 
using Equation (1). Table 1 lists a = 1.3 DELX as the 
maximum likelihood estimate. An investigator would hope 
that if the true threshold variance were equal to 1.00, 
then DELX would be close to .77* If so, then a would 
appear unbiased and qualify as a useful input to follow-up 
tests. However, the average length of DELX from 1000 Si 
sequences generated from N(0, 1) sensitivity data equalled 
.49. With an S 2 or S 3 sequence, an investigator would 
hope that DELX would be close to .33, since a = 3 DELX. 

Again a large sample (675) was investigated with an average 
DELX of .21. Thus the maximum likelihood principal 
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(Equation (1)) yielded an average estimate of a that was 
substantially smaller than the true standard deviation. 
Figures 10 through 12 are histograms of the frequency 
distributions of DELX for Si, S 2 and Sa . Table 2 lists 
the same data. 

A summary of the associated a’s from the simulated 
sequences was as follows: 





a < a/2 


a/2 


< a < 2a 


a > 2a 


S 1 


45.2# 




50.6% 


H.2% 


S 2 


51. W 




m.1% 


0.9% 


S 3 


51.15S 




1 A% 


1.2% 


The 


estimates of a were improved slightly 


by using 


an extension 


of three trials 


of S 


1 . The results 


were as 


follows . 


a < a/2 


a/2 


< a < 2a 


Q > 

V 

rv) 

Q 




W 




W 


2 % 



Figure 13 is a histogram of the latter distribution. 
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RANGE OF DELX 



Si 



S 2 



S 3 



.00 


to 


.02 


29 






26 


”1 




17 


.02 


to 


.04 


24 






32 






28 


.Oil 


to 


.06 


26 






30 


1 




24 


.06 


to 


.08 


34 






14 


Q> 

A 

Q 

ro 


18 


.08 


to 


.10 


7 






5 






0 


.10 


to 


.12 


53 






48 






43 


.12 


to 


.14 


57 






25 






35 


.14 


to 


.16 


0 






0 






0 


.16 


to 


.18 


0 


o<a/2 


0 


— 




0 


.18 


to 


.20 


6 






8 






1 


.20 


to 


.22 


35 






17 






22 


.22 


to 


.24 


84 






52 






36 


. 2 -4 


to 


.26 


68 






29 






30 


.26 


to 


.28 


17 






6 






11 


.28 


to 


.30 


0 






0 






0 


.30 


to 


.32 


0 






0 






0 


• 32 


to 


.34 


0 






0 






0 


.34 


to 


.36 


0 






0 






0 


• 36 


to 


.38 


1 






0 






0 


.38 


to 


.40 


13 






1 


o/2<a<2a 


2 


.40 


to 


.42 


24 






3 






5 


.42 


to 


.44 


46 






13 






8 


.44 


to 


.46 


46 






15 






14 


.46 


to 


.48 


63 






8 






9 


.48 


to 


.50 


44 






6 






6 


.50 


to 


.52 


32 






4 






4 


• 52 


to 


.54 


23 






5 






8 


.54 


to 


.56 


12 






0 






0 


• 56 


to 


.58 


1 






0 






0 


.58 


to 


.60 


0 






0 






0 


• 

• 


• 

• 


• 


• 

• 


a/2<a<2a 


• 

« 




< 


• 

• 


.72 


to 


.74 


• 

0 






• 

0 






• 

0 


.74 


to 


.76 


1 






0 






0 


• 76 


to 


.78 


0 






0 






0 


.78 


to 


.80 


6 






0 






0 


.80 


to 


.82 


8 






0 






0 


.82 


to 


.84 


11 






0 






0 


.84 


to 


.86 


13 






0 






0 


.86 


to 


.88 


23 






0 






0 


.88 


to 


.90 


19 






1 






1 


.90 


to 


.92 


18 






1 






1 


.92 


to 


.94 


21 






0 






0 


.94 


to 


.96 


23 






0 






0 



TABLE 2 

FREQUENCY DISTRIBUTION OF DELX 
STIMULUS ~ N( 0 , 1 ) 



1 

o<o / 2 

J 

a/2<a<2o 
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RANGE OF DELX 



S 3 



.96 


to 


.98 


20 


0 


0 


.98 


to 


1.00 


11 


0 


0 


1.00 


to 


1.02 


10 


0 


0 


1.02 


to 


1.04 


10 


• 


• 


1.04 


to 


1.06 


6 


• 


• 


1.06 


to 


1.08 


6 


• 


• 


1.08 


to 


1.10 


4 






1.10 


to 


1.12 


3 






1.12 


to 


1.14 


0 






1.14 

• 


to 

• 


1.16 

• 


0 

• 






• 

1.56 


• 

• 

to 


• 

1.58 


• 

0 






1.58 


to 


1.60 


1 






1.60 


to 


1.62 


1 


o>2a 


1.62 


to 


1.64 


2 






1.64 


to 


1.66 


0 






1.66 


to 


1.68 


0 






1.68 


to 


1.70 


2 






1.70 


to 


1.72 


1 






1.72 


to 


1.7** 


1 






1.7** 


to 


1.76 


1 






1.76 


to 


1.78 


3 






1.78 


to 


1.80 


0 






1.80 


to 


1.82 


0 






1.82 


to 


1.84 


5 a>2o 




1.84 


to 


1.86 


2 






1.86 


to 


1.88 


2 






1.88 


to 


1.90 


2 






1.90 


to 


1.92 


5 






1.92 


to 


1.94 


3 






1.9** 


to 


1.96 


1 






1.96 


to 


1.98 


2 






1.98 


to 


2.00 


2 






2.00 


to 


2.02 


2 






2.02 


to 


2.04 


1 






2.04 


to 


2.06 


2 


• 


• 


2.06 


to 


2.08 


0 


• 


• 


2.08 


to 


2.10 


1 


• 


• 


2.10 


to 


2.12 


0 


0 


0 


2.12 


- — 




0 _ 


0 _ 


0 






TOTAL 


1000 


350 


325 



a>2o 



TABLE 2 (Continued) 
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V. CONCLUSIONS AND RECOMMENDATIONS 



Estimates of y using MBS are suitable for follow-up 
tests while those of a are not. Additional trials will 
be necessary to refine the estimate of a. 

Termination of MBS with a small number of trials cannot 
be guaranteed. It is not possible for an investigator to 
determine beforehand the number of trials necessary to 
complete a terminal situation (Si, S2 , or S3). The average 
number generated in this investigation was eleven with 
extreme values of five and twenty-eight. It appears that 
this number is a function of the location of endpoints, 
the variance of the threshold distribution, and the sequen- 
tial outcomes of the experiment. 

The extension procedure in MBS should be investigated 
so that levels yielding the most useful data can be tested. 
For example, the middle level tested in an Si extension 
doesn't offer data to confirm or alter previous estimates 
of y and a whereas other levels may do so significantly. 

As an alternative approach to obtain practical estimates 
of c, a Monte Carlo method is suggested in place of the 
maximum likelihood principle. Note that sample values of 
DELX are "bunched" together in distinct intervals separated 
by intervals in which no value of DELX occurs (see Table 2 ). 
If an investigator could determine the interval to which 
the DELX from his experiment belongs, then he should be 
able to estimate o better than before. 
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COMPUTER SEARCH FOR MLES 



THE FOLLOWING PROGRAM FINDS THE MAXIMUM LIKELIHOOD 
ESTIMATES OF THE MEAN AND STANDARD DEVIATION OF ANY 
SEQUENCE OR EXTENSION. 



P SI ( 1 ) IS THE INITIAL ESTIMATE OF THE MEAN. 

PSI ( 2) IS THE INITIAL ESTIMATE OF THE STD. DEVIATION. 

DELLC IS THE DESIRED ACCURACY OF THE MLES. 

MAXEV IS THE MAXIMUM NUMBER OF ITERATIONS ALLOWED 
FOR ANY SEARCH. 

DELCAP IS THE INITIAL STEP SIZE. 

RHO IS THE REDUCTION FACTOR OF THE STEP SIZE ONCE A 
LOCAL MAXIMUM HAS BEEN FOUND. . 



883 PSI ( 1 ) =T( NL) +DELX 
885 PSI ( 2 ) =2. 5*DEL X 
DELCAP=DELX 
DELLC=. 05*DELX 
RHO=. 12 5 
MAXEV=75 
SLC ( 1 ) =DELCAP 
SLC( 2 ) =DELCAP 
SPSI=S( PSI (1) , PSI ( 2) ) 

I VAL=1 

801 SS=SPSI 

PHI ( 1 )=PSI ( 1 ) 

PHI ( 2 ) =PS I ( 2 ) 

I BK=- 1 
GO TO 840 

813 IF( SS.LT.SPSI ) GO TO 803 

802 I F ( IVAL.GE. MAXEV) GO TO 852 
DO 820 K = l,2 

IF C SLC ( K) I 821 » 852 t 822 

821 IF(PHI ( K) .LE.PSI (K) ) SLC ( K ) =-SLC ( K ) 
GO TO 823 

822 IF(PHI(K).GE.PSI(K) ) SLC ( K ) =-S LC ( K ) 

823 THET = PSI ( K) 

PSI ( K) = PHI (K ) 

PHI ( K) = 2.*PH I ( K) - THET 
820 CONTINUE 

I F ( PH I ( 2 ) ) 815,815,816 

815 PHI (2)=PSI (2)/2. 

816 SPSI = SS 

SPHI = S( PHI ( 1 ) , PH I ( 2 ) ) 

SS=SPH I 
IVAL=IVAL + 1 
I BK=1 

840 DO 841 K=1 , 2 
THET = PHI ( K ) 

SLC I = SLC(K) 

PHI(K) = THET + SLC I 
I F ( K-2 ) 824,825,825 
825 IFIPHI ( 2) ) 829,829,824 
829 PHI ( 2) =THET/2. 

824 SPHI=S( PHI ( 1 ) , PHI ( 2) ) 

I VAL = I VAL + 1 

I F ( SPHI .GE. SS) GO TO 842 
PHI ( K ) = THET - SLCI 
I F ( K-2 ) 818,819,819 

819 I F ( PH I ( 2 ) ) 817,817,818 
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817 PH I ( 2 ) = THET/2. 

818 SPHI =S ( PHI (1) , PH I ( 2 ) ) 

I VAL = I VAL + 1 

I F ( SPHI.LT.SS) GO TO 844 
SLC(K) =-SLC I 
842 SS=SPH I 
GO TO 841 
844 PHI (K) = THET 
841 CONTINUE 

IF(IBK) 813,813,828 
828 IF(SS.LT.SPSI) GO TO 801 
DO 826 K=1 , 2 

IF(ABS(PHI (K)-PSI(K) ) . GT . 0. 5* ABS ( S LC ( K ) ) ) GO TO 802 
826 CONTINUE 

DELF=DELCAP 

803 IFIDELCAP.LT. DELLC) GO TO 852 
DELCAP=RHO*DELCAP 
SLC( 1)=RH0*SLC( 1 ) 

SLC(2)=RH0*SLC(2) 

GO TO 801 



STATEMENT 852 IS THE EXIT FROM THE ROUTINE. 

PHI ( 1 ) IS THE FINAL ESTIMATE OF THE MEAN. 

PHI ( 2 ) IS THE FINAL ESTIMATE OF THE STD. DEVIATION. 
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